V «. i' - v 


™mL\ »* 

-* 

?!-'•*•' .3 : '■ 


v-M * v r;vm ? 

.•'. ■ l 1 !i v.;; y • ;;>. 


•■ ; cvv'iir .* ,•?> 

JM ' . 


i «, 


V. :.. '■; ■ V 

*-' . '-f * ! \-\) ,: • ' . Vj '-Vi 

; ■ i ■ .' v ' * 

• ^ }■ • - .:■', . ? J > 


. ■ v 

’ * 

. r. 
v ■• , 









'\ •f; /* 


’ I 1 ./ 






;#/ V 

: V*w.V • r • , ■ 


I' 1 -® 


|l| '' , - . . ^ Whi rtir thp - \ * 

^*- s /:■ • • -A ••/ . .;■■ , f ’■ .£■•' •’.-A •-•'/ rzv?w ’* 'vPWw w jv- -r. v- •■ <*- ./ •.;*./¥•: i;:-v*v :•-•/ 

E. . • ■ ■ ■ : « FW"ACO^|C JN PUCTS-;. ; ;, . ■ , . , 

- , - . ^ - r^ ' \ ^ y I'.’i? ; vO ^ ^ ^ ^ ^ ^ ^ ^ ^ 1 1 


‘ - •■/// 

- • - • . - . ‘ ■ 

.. - . .. ' .V;w i. : ... I , .. 

fegs i : .>v . - •;.. V v{-r-:.v ;.< '. . Mi ■■■I-..-: 

V^Vi ,v *V- -’■ 'V L— . 



W&A 


,$#v; 

■i.v-v . 





,.;h J -V ■ -• '*■'" : ■ ' - ... ■ .• V > ••.>••.•/ ■"’•■•. 

W. ' ; • : - ( - ■ 1 y^s - *. •_ 

|| i- | . \ r 


i if; • . 


k.i/ / • ■■• w . . ‘ 


't: 





(NASA-CE-179845) I0NEAMZN1AL INVESTIGATIONS 
OF TBE FINITE ELEMENT SCLDTIGNS FOB ACOUSTIC 
FEOBAGATICN IN DUCTS Final Eepoct (Missouri 
Ociv. ) 86 p CSCL 20A 


N87-1Q751 


Unclas 

G3/71 44333 


Iff. 


llfc 


: ; - 


University of Missouri-Rolla 






■ '$t; a :: y V sAyffy- ■ 'y ■ v - p '^ 




■■ ■ '> i-^r.nvjs;.; -w’.-V 


i I •' V- 

•<;•: •'■•' ./: .7 ;T ^ >;‘T 

,.•'•• !4V. .a'- :»• :v,- 


FUNDAMENTAL INVESTIGATIONS OF THE FINITE ELEMENT SOLUTIONS 
FOR ACOUSTIC PROPAGATION IN DUCTS 


N. J. Walkington*, T. M. Wickst and W. Eversman* 
University of Miss our i-Ro 11a 
Rolla, Missouri 65401 


FINAL REPORT 
NASA LANGLEY NAG- 1-1 98 


*Department of Mechanical Engineering 
t Department of Engineering Mechanics 



ABSTRACT 


The question of convergence of three finite element algorithms for 
the modelling of acoustic transmission in ducts carrying a nonuniform 
mean flow is addressed. The details of each algorithm are stated and 
example calculations in uniform and nonuniform ducts are made and assessed 
for accuracy and convergence. 

The algorithm based on the assumption of irrotationality is found to 
be highly convergent. This algorithm is the one used in current turbo-fan 
inlet acoustic radiation codes. A theoretical analysis indicating con- 
vergence is supported by example calculations. 

Two additional algorithms which do not require irrotationality are 
found to be less convergent, and perhaps not convergent at all for certain 
severely sheared velocity profiles. No theoretical convergence criteria 
can presently be established for these algorithms and convergence difficul- 
ties are shown here by example. Included in this class of algorithms is 
the duct analysis program ADAM which is known to display apparently non- 
convergent solutions in certain cases. 
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I. INTRODUCTION 


It has become apparent that the several finite element approximations 
used in the modelling of acoustic propagation in nonuniform ducts do not 
possess consistent convergence characteristics. In particular, experience 
with approximations based on an assumption of irrotational mean flow and 
acoustic perturbation indicates that, at least over the fairly broad 
parametric range so far investigated, convergence characteristics are 
good. On the other hand, when the irrotational assumption is not imposed, 
necessitating a finite element approximation based on the primitive form 
of the acoustic equations (in terms of pressure and two or three velocity 
components), serious convergence problems can occur, particularly when 
the mean flow has steep gradients, representative of typical boundary 
layers. It is the purpose of this study to attempt to isolate with 
simple examples the nature of the difference in convergence characteristics 
of the finite element algorithms. 

In modelling acoustic propagation in ducts the following assumptions 
are often made: 

1. The fluid is an inviscid ideal gas. 

2. The linearized equations of motion accurately describe 
the propagation of the acoustic disturbances. 

3. The motion of the fluid is irrotational. 

For the aircraft engine inlet problem assumptions 1 and 2 imply assumption 
3 since, for this problem, the fluid is initially irrotational. By the 
Helmholtz Vorticity Theorem, the fluid remains irrotational. However, 
it is known that near surfaces assumption 1 breaks down, since the 
effects of the fluid viscosity are then significant. In order to com- 
pensate for this deviation, without giving up too many of the simplifica- 
tions afforded by assumptions 1-3 the following could be assumed: 
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4. The motion of s mall disturbances is well described by the 
linearized inviscid equations with a viscous mean flow 
substituted in place of the inviscid mean flow. In other 
words, the coefficients of the linear equations are altered 
to represent a viscous mean flow. 

Assumption 4 has no physical, mathematical or intuitive justification, 
and for this reason it is generally argued that it is not appropriate. 
However, we do make some calculations with it here to demonstrate con- 
vergence. 

Several numerical algorithms have been developed to model the duct 
acoustics problem, [1-14] , probably the most versatile of which are the 
finite element methods (FEM's). However, it has been reported [7] that 
in certain instances some algorithms do not converge. In this document, 
we consider three finite element algorithms used to model duct acoustics 
problems; one based on a velocity potential formulation and the other 
two based upon the primitive equations. It is shown that the algorithm 
using the velocity potential will converge for a very broad class of mean 
flows and geometries. The other two algorithms are not ammenable to such 
analysis, so numerical experiments were conducted in order to estimate 
what problems may be solved with these algorithms and to compare the 
algorithms . 

Fix, Nicolaides, and Gunzburger [11-13] consider a least squares 
algorithm to solve the acoustics problem using the primitive equations. 

They prove that for the no flow case their algorithm will converge optimally 
when a special class of finite element grids is used. A similar least 
squares algorithm was proposed by Astley and Eversman [3]. However, they 
reported that better results were obtained with a classical Galerkin type 
approximation . 



II. ACOUSTIC EQUATIONS AND BOUNDARY CONDITIONS 


In this section we consider the acoustic equations and boundary 
conditions for a duct, ft. 

Assumptions 1 and 2, or 4 lead to the following equations for 
the acoustic velocity and pressure in ft. 


p + V.(pu + ^ p U) - 0 
c c 

iwu + (U.V)u + (u.V)U + V(^) * 0 

-V ~ -s, "O' P 

where 


( 2 . 1 ) 


p,u are the acoustic pressure and velocity 

p,c,U are the mean flow density, speed of sound and velocity, respectively, 
a) is the frequency. 

If assumption 1 is used the mean flow (p,c,U) should satisfy 


V.(pU) - 0 

(U.V)U + - V (tt) - 0 
P ~ 

TT - < P \ c 2 

* P 


( 2 . 2 ) 


where 


Y and < are constants and y ■ c - c . 

1 p V 

When assumption 4 is used the mean flow should satisfy the viscous 
compressible equations. However, these are difficult to solve so often a 
boundary layer correction to an inviscid flow is used as a further approxi- 
mation/assumption. 



If the irrotational assumption is used equations (2.1) reduce to 


7 + !-(p^+\p?) =o 

■jj p + iu)<p + U.V<J) * 0 

where 

4> is the acoustic velocity potential. 

If assumption 1 is used the mean flow will satisfy, 

7(p7|) - 0 

1 

p = [p Y_1 - V$.V$] Y_1 

where 


(2.3) 


(2.4) 


$ is the mean flow velocity potential 

A 

p is the (constant) stagnation density. 

A 

When assumption 4 is used (p,c,U) represents the rotational mean 

flow. 

BOUNDARY CONDITIONS 

Since we are interested in the duct problem we consider the domains, Q, 
to be of a specific form. We assume (see figure 1), 

an = r l u r 2 u r 3 u r 4 (2.5) 
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where- 


T l = { (0 ,y) : 0 < y < d^} 

* left hand end of the duct 
r 2 = {(£,y): 0 < y < d 2 } 

■ right hand end of the duct 
r 3 « lower duct wall 
r 4 - upper duct wall 

d^, d 2 are the duct widths at each end and l is the duct length. 

We assume that the duct walls are "hard". This requires that 


or 


u.n ■ 0 and U.n ■ 0 


|i - 0 and 
9n 9n 


on 


r 3 u r 4 


( 2 . 6 ) 


where n is the unit outward normal to the duct. 

It is interesting to note that when assumption 4 is used the mean 
flow will satisfy U * 0 on the walls while the tangential acoustic velocity 
may or may not vanish. 

It will be assumed that the duct terminates with semi-infinite 
uniform ducts which carry a uniform plug flow. In this situation, the 
solution in the duct, ft, may be matched on to an analytical solution in the 
semi- inf inite ducts. This solution, for the irrotational problem, is given by 

OO 

<J> = £ [a* exp(A + x) + a” exp (A x) ] cos(^j y) (2.7) 

n=0 n n n n d 


6 




) 
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where 


d is the duct width 

{a + }, {a - } are arbitrary constants 
n u 



n 



(ikM + 6 ) 
n 


M » | U j /c » Mach, number of the mean flow 
k » ^ “ spatial wave number 

6 - i[k 2 - (l-M 2 )^) 2 ]* 5 

n a 


When rotational disturbances are permitted the uniform duct solution is 


Z 

n"0 


Cos (££ y) 


Gos(j^ y) 


Sin(|^y) 


[A ] 
n 


a + exp(X + x) 
n r u 

a" exp (X x) 
ti n 

0 ,,0 . 

L*n expa n 


( 2 . 8 ) 


where 


tv 


- X" 


pc (ik + X^M) 


1 


(?) 

d t * , . “h 


pc(ik + tt M) 
n 


- X 


pc (ik + X M) 
n 


<f> 


pc(ik + X M) 
n 


M / HTT . 
k ^ k' 


{a + }, {a - } and {a°} are amplitude coefficients 
n n n 


X~ are as in equation (2.7) 

X° --$i 
n M 

The last column of the coefficient matrix [ A^] introduces the rotational part 

of the disturbance. Similar representations are used at x * l, with [ A^] 
replaced by [1^] and the a n replaced by b n . 




The quantities {a + } and {a } represent wave amplitudes of waves moving 
n n 

in the positive and negative directions respectively. The {a||} are the 
amplitudes of hydrodynamic waves that move with the mean flow and do not 
exist if the flow is irrotational. At the duct terminations we require 
the solution within the duct to match the above analytical solution. 

In particular we require 

OO 

(PV<|> + , 4up). n | » p. Z <5* (a* ” a”) Cos(~ y) (2.9) 

c ~ (0 ,y) 1 n-0 n n n d l 


OO 

W + ij 0 p). n| (i>y) - -02^ 6*(b+ - tT) Cos(SI y) 


or, when potential flow is not enforced 


(pu + p U). n | » i [. 


<5 


(0,y) n-0 cCik+X^M) 311 c(ik + A~M) ~ n 

n n 


n 


PM .mr. 0 , ,mr . 

+ k <3^ a n Co3( d7 y) 


( 2 . 10 ) 


(pu + p U) .nj 
- c ~ ~ (*,y) 


Z [ 

n-0 c (ik + A^M) 


b + 
+..v n 


n , - , PM yttir „ , 0 , „ ,mr ^ 

b + — (j-)b l r Cos(-r- y) 

c(ik + A M) k d 2 V 2 d 2 


Where in equations (2.10) the subscripts 1*1 and ?2 indicate that the quan- 
tities appearing within the brackets are to be evaluated on and T 2 respectively. 
We specify the waves incident upon the duct, fl. This corresponds 

to specifying {a"*"}, (b } and {b^} if the mean flow Mach number is negative 
n n n 

(specify {a|j} if M > 0). 

If the duct terminations are not in regions of uniform plug flow, 
no analytical solutions corresponding to equations (2.7) and (2.8) exist. 
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Astley and Eversman [14] proposed a method for generalising these modal 
boundary conditions. However some technical problems arise with this 
approach so only plug flows are considered here. 
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III. WEAK PROBLEMS 


In this section we discuss several weak formulations used for the 
modelling of the duct problem . We first discuss a formulation using 
an acoustic velocity potential. The other two weak formulations, based upon 
the primative form of the equations, are those pursued by Astley and 
Eversman [2-3] and Abrahamson [6-7] respectively. 

1. Velocity Potential Formulation 

In order to incorporate the modal boundary conditions of equation 
(2.9) into the problem it is necessary to define a space of sequences, S, 
which are square summable in a weighted sense. Let 


S « {a: a 2 + Z |a | 2 n < °°} 

U n-1 n 


(3.1) 


The modal amplitudes at the left and right hand ends of the duct will be 
denoted by (a + ,a ) and (b + ,b“) respectively and will be required to lie 
in S x S (i.e. a + eS and a“£S, etc.). The sequences of eigenvalues of 
the ends of the duct are denoted by (X + ,X ). In formulating the weak 
problem we will need certain spaces. First recall that L^(ft) defines the 
linear space of equivalence classes of measurable functions g such that 

/ |g| < 00 and that H (ft) defines the linear space of equivalence classes 

2 2 
of functions in L (ft) whose first distributed derivatives are also in L (ft) . 

The weak problem may now be stated. Find (<p,p,a ,b ) £ H (ft) x L (ft) x S x S 

such that 


(i) { t P t|) a (pV4> + p U).V i|>} 

J C C “ - 

ft C c 


“ P 1 2 ^ a n “ O C0S ^ ^ 

1 n n n a. 


n=0 


l l 


+ p 2 2 o> ■ b ) 6 c °s(t- y) t 

2 n-0 n n n l d 2 


(3.2a) 


(ii) 



U.V<J> + iuxp }q 


0 


0.2b) 
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(ill) Z n{|-(a + + a j 
n-1 d l n n 


* C03(SIy)}c. n+ {i- (aj + V 


4>}cx rt “ 0 


(3.2c) 


(iv) Z n{f-(b+ + bj 
n-1 d 2 n n 


♦ Cos(f^y)}p n + {i- (-J + aJ) 


<j)}6 n - 0 y (Ij»,q,a,6) 6 H 1 (ft) x L 2 (G) x S x S 


(3. 2d) 


In the statement of this weak problem (p,c,U) is the mean flow, a + and b” 

are presumed to be given and p^, dj^ and P 2 » are the values of the mean 

flow density and duct widths at the left and right hand ends of the duct 

1 2 

respectively. The quantities 6 and 6 are the radicals defined in 

n n 

equation (2.7), evaluated at the left and right hand ends of the duct 
respectively. 

It is shown in Appendix I that this weak problem is well posed under a 
physically reasonable hypothesis on the mean flow. 

We have chosen to solve for the pressure explicitly. As an alter- 
native the velocity potential could be calculated initially by eliminating 
the pressure and then the pressures calculated subsequently. 

2. Astley Eversman Formulation 

The modeling of the modal boundary conditions (equations (2.9-2.10) 
is more involved when the primative equations are utilized. Astley and 
Eversman proposed the following ideas to accomplish this. 
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For each (p,u) 6 H d (ft) x H^(fl)^ define six sequences a + , a , a^. 


b + , b , b^ where 


a 

n 

_ i 

• 

,mr . “1 

cos (-p y) 


p 

a 

n 

i * 2 

- [A 3 j- j 

n d t j 

1 

1 

,nlT x 
cos (-3— y) 

d l 


u l 

0 

a 

u n -J 



. cos (t~ y)_ 

a l 


- U 2 - 


( 3 . 3 a) 


rb + i 

n 

-1 


r /HTT . 
cos (7- y) 

a 2 


P 

b~ 

n 

- [B ] j 

n d 2 j 

1 

2 

.nir v 
cos (j - y) 

a 2 


U 1 

b° 

■- n J 


4 

.mr . 
cos (-5— y) 

l a 2 


- u 2 - 


( 3 . 3 b) 


for n=l, 2,3 with ajj * bg * 0 


r c 


L a 0 J 




L»OJ 


- a; 


i 

- a: 


pc(ik + AqM) pc(ik + A“M) 


-1 


- a; 


1 

- a: 


i 


pc (ik + AqM) pc(ik+A Q M) 


ti [;] 


d. 


r 2 ‘2 


([:,] 


( 3 . 4 a) 


( 3 . 4 b) 


We now fix some notation so that the space of solutions may be defined. 
For w = (p,u) e H 1 (fl) x H 1 (fi) 2 define 


A + : H 1 (fl) x H 1 (ft) 2 + l 2 


by 


A + (w) * a + 

where a + is defined in equation ( 3 . 3 ). Similarly we difine A - , A 0 , B + , 
B , and B^. 
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Let a + , b~, be specified. The function spaces of Astley 
and Eversman are: 

U « {w 6 H 1 (£2) x H 1 (£2) 2 : A + (w) ■ a + ,B”(w) ■ b”, 

B°(w) - b 0 } 

Uq=* {w 6 (£2) x H 1 ^) 2 : A + (w) * B~(w) « B^(w) - 0} 

with the weak problem becoming: 


Find (p , u) 6 U such that 


(i) 


{ pq + (pu + -^Dp). Vq} 


a 


.1 


- z {- 


a - 


n«0 c(ik + X + M) n 
n 


- P 1 M „ n 

2 a' + -f- (fl) a 0 } 

c(ik + jfM) “ k d l “ r, 

n i 


Cos (~ y) q 
a l 


+ Z {- 


B— b + 

+wx n 


n=0 c(ik + M) 
n 


b n + <?> b n } 

c(ik + ~M) K d 2 T_ 

n z 


Cos(B^ y) q 
a 2 


- 0 


(ii) j (iwu + (U.V)u + (u*7)U + V(*)}. v - 0 

n 

V (q, v) e U 

Here it is understood that a + , b and b^ are data with a = A (p,u), 
b * B + (p,u), and « A^(p,u). 


(3.6) 

(3.7) 

(3.8a) 


(3.8b) 
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3. Abrahamson Weak Problem 


Abrahams on [6] used Che following weak problem for his ADAM 
computer program. Lee 

V" { (p » u) 6 U: u.n| r * u.nL * 0} (3.9) 

~ ~ A 3 - - *4 

v 0 -vnu 0 . 

then find (p,u) e V" such that 

{ij p + V. (0 u + U p)}q - 0 (3.10a) 

a' ' ' c 

(iu>u + (U.V)u + (u.V)U + V(£)}.v - 0 V(q,v)eT^. (3.10b) 

J *• ~ ~ ~ ~ ~ ~ ~ P ~ u 

a 
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IV. FINITE ELEMENT APPROXIMATION OF THE WEAK PROBLEMS 


1 . Velocity Potential Formulation 

1 2 

The function spaces H (ft) and L (ft) from which <j> and p are chosen 
may be approximated by the usual f inite element subspaces , (ft) and 1M^ (ft) . 
The superscript is used to indicate the degree of inter-element continuity 
(0 continuous, -1 discontinuous) and the subscript indicates the degree of 
polynomial used in the construction of each element. To approximate the 
sequences the following finite dimensional subspaces are utilized. 


S„ * {a 6 S: a * 0 v n > N} 
N n 


(4.1) 


The approximate weak problem is then; Find (<p,p,a ,b + ) e M^(ft) xM 
x Sjj such that equations (3.2a-d) hold V 0J>,q,a,S) eiM^(ft) x 
]M _1 (ft) x S N x Sjj. 


-1 

l 


(ft) x 


Under certain technical assumptions it is shown in Appendix 1, that 
if > 

(i) a > k-1 

(ii) N > zr where h is the maximum diameter of the elements used 

— n 

in the finite element triangulation, 

then. 


I* - * 0 ii - c i hk 


(4.2) 


Ip “Polo — c 2 hk ^ c l* c 2 constants ) 

where (4 > q»Pq) is the solution and (4>,p) is the approximate solution. 

In other words the finite element scheme converges with an optimal rate 
if suitable care is given to a required compatibility between the approximation 
for the velocity potential and the pressure and boundary coefficients. 
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2 . Astley-Eversman Algorithm 


It is not immediately obvious how the class of functions U 
discussed in section III. 2 may be approximated. Astley and Eversman 
[3] proposed the following algorithm to generate (non-conforming) 


approximations . 

Let (ft) 6 H^(£2) x be the usual finite element 

subspaces as described in section IV. 1 above. On the trace of any 
element (u,p) 6 M^(&) x may be written as 





correspond to the function values of 


ta 


at node i on 


are the interpolation functions corresponding to node i on The 

nodal values of (p,u) are then constrained to satisfy, 



(4.4) 
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for all j » 1,2, ... M. N represents the number of modal amplitudes 
being used in the approximation, and [A^] is the matrix defined in 
equation (2.8). Equation (4.4) gives a relationship between the modal 
amplitudes and the nodal values on the boundary for the approximation. 
This relation may be written 


Ho ’ [T o ] 6 


<4. 5a) 


where Uq corresponds to the nodal values of (p,u) on and A the vector 
of nodal amplitudes (a*, a^, a^, n * 1,2,...N). An analogous expression 
for the nodal values on may be written as 


Hu - IV ?• 


(4.5b) 


If a typical element of M^(ft) x M^(£2)^ is described by 


"P - 


rnoi 

U 1 

- m 

u 

lu 2 J 


.Hz. 


(4.6) 


where Uq and are the nodal values on F^ and and U contains the 
other nodal values a typical element of the space used to approximate 
U would described by 


P 

U 1 

- [*i 

i 

M ? 
O 

e- f 

i 


"a' 

u 



i 

i H 

1 


_B_ 


(4.7) 


where I is an identity matrix and [Tq], [T^J are as in equation (4.5). The 
constraints on a + , b , etc. can now easily be implemented using penalisation. 
Astley and Eversman [7] describe the implementation of this algorithm in 
detail. 
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V. RESULTS 


1. Uniform Duct Results 

The analytical solution to the acoustic problem for a uniform duct 

carrying a uniform mean flow was given in equations (2. 7-2. 8). This provides 

a test case which can be used to estimate how well the finite element approx- 

2 

imations do approximate the solution. Throughout this section we use the L 
norm of the error in pressure, defined by, 

2 

p - p h i 

ft 

where 


P - P 


h 1 0 


p * analytical solution 
p^ 3 finite element solution, 

to estimate the error in the solution. Standard results [15] show that a 

finite element approximation using polynomials of degree k would have an 

2 

optimal rate of convergence in the L norm of order k+1, i.e. 



0 i c h 


k+1 


where 


h * diameter of largest element, 

if the solution sought was sufficiently smooth. Fix, Nicolaides, and 
Gunzburger [11] obtained optimal convergence for the no flow case using 
their least squares scheme with linear triangles (k=l) on a "criss cross" 
grid (i.e. a grid satisfying the grid decomposition property [10]). 


I 
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Potential Flow Code: 


We consider the potential flow code separately from the other 
two algorithms for two reasons. First, we have the analytical results 
in the Appendix that indicate that the finite element scheme should 
converge optimally when quadratic elements are used for the velocity 
potential and discontinuous linear elements are used for the pressure. 
Secondly, the results indicate that extremely ’rough' mean flows are 
admissible since no mean flow derivatives enter this problem. 

When the finite element approximation described in equations (3.2) 
was used to solve for the analytical solution given in equation (2.7) 
optimal convergence of the pressure was observed in all cases tried. 

This was expected since the analytical solution has infinitely many 
derivatives. A more interesting example may be constructed by combining 
two of the solutions given in equation (2.7) with different mean flows. 
(See Figure 2) This results in a problem with a discontinuity in the 
mean flow and the acoustic pressure, the velocity potential being 
continuous. The compatibility condition at the interface may be 
written as 

£(pu + Up) • ej (0.5) » 0 
c 

We consider the results using the following as data: 

co = 1.0, a^ * 1.0 a + = 0, n=l,2,..., £=1.0, d = 0.5 

u n 

c = 1.0, bg = 0 n = 0,1,2,.... 

v 1.0 if x < 0.5 TT _ -0.25 e v x < 0.5 

P(x,y) = U = ~x 

0.5 if x ^ 0.5 ~ -0.5 e x x ^ 0.5 
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The choice of a plane wave incident at x*0 (a+ ■ 1) is presented so 
that the jump in pressure at the interface can easily be plotted. 
Calculations were made with non plane waves and the results were 
essentially the same. 

Figure 3 shows the real part of the pressure plotted near x=0.5 
that was calculated using a mesh with 10 uniformly spaced elements across 
the duct and 20 elements along the duct (a 10 x 20 mesh). This solution 
almost exactly overlays the analytical solution. It can be seen that 
the jump is captured without any of the characteristic wiggles associated 
with ill-posed problems. It is interesting to note that while the 
pressures away from the discontinuity could be discontinuous, the solution 
has almost no jumps at the other element interfaces. The 10 x 20 mesh 
used for figure 3 has element boundaries along the interface x=0.5, so 
the jump can be captured exactly. If a 10 x 21 element mesh is used 
the discontinuity can no longer be captured exactly since the jump in 
pressure does not coincide with element boundaries. The solution for 
the real part of the pressure near x*0.5 obtained with a 10 x 21 mesh 
is given in figure 4. It can be seen that in this situation the pressure 
in the element centered at x=0.5 simply stradles the jump and otherwise 
the pressure is still very accurate. It is possible to use elements 
with continuous pressures to solve this problem. This was done for both 
the 10 x 20 and 10 x 21 element meshes and the pressures are plotted in 
figure 5 and 6. The error associated with the discontinuity appears to 

spread over at least three elements for each of these solutions. 

2 

A plot of the L error in the pressure versus mesh size (h) is given 
in figure 7. When discontinuous elements were used that had boundaries 
on x=0.5 the convergence rate is optimal (slope = 2). The other two lines 
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with slope * 1 are the convergence curves for the case where the inter- 
face x*0.5 always bisectes an element so that the jump couldn't be 
calculated exactly and the case where the pressures were continuous. 

Since this problem does not have a smooth solution the approximation 
results given in the Appendix do not apply and this explains why even the 
discontinous pressures may give sub- optimal convergence rates. 


ASTLEY AND ABRAHAMSON CODES: 


The weak problems presented in equations (3.8) and (3.10) both 
involve derivatives of the mean flow so discontinuous mean flows, like 
the one discussed for the potential flow code, can not be considered. 

In order to investigate the behavior of the two algorithms, a uniform 
duct case was considered; the solution being given in equation (2.8). 
Since the solutions were all continuous it is possible to calculate the 
error of the solution defined by: 


2 


p - p 



9 

t tip - P h | 2 + |V(P - P h )| 2 } 

h 


where 


p » analytical solution 
p^* finite element solution 


Standard results show that if the solution sought is smooth, and 
polynomials of degree k are used to construct the finite element mesh 
approximating the solution, the optimal rate of convergence is k, i.e. 

I P - P h I 1 c h k 
n 1 

2 1 

Error calculations were made using both the L and H norms. 


22 



The following data are used for the uniform duct problems considered. 

co » 1.0, l ■ 1.0, d = 0.5 

U ■ -0.5 e p » 1.0, c » 1.0 
~x * 

a„ * b n =■ 0, a + =b^=0 n * 2,3,... b =0 n = 0,1,... 

u u n n n 

Two cases are discussed, the first has a wave of unit amplitude incident 
at x**0 corresponding to the first transverse mode (a* = 1.0), and the second 
case considers a wave of unit amplitude incident at x=l .0 corresponding 
to the first hydrodynamic mode (b^ * 1.0). This latter case has no 
correspondence in the irrotational model so the primitive form of the 
equations must be used to solve for this solution. For each of the two 
algorithms eight noded isoparametric elements were used for both the 

pressure and velocity. Optimal rates of convergence would then be 2 in 

1 2 
the H norm and 3 in the L norm. 

Figures 8 and 9 are plots of the and errors in velocity and 

pressure respectively, verses mesh size for the first problem. For this 

problem both algorithms gave (approximately) the same convergence rates, 

however these rates are not optimal. The rates were one and two for 
2 

and L errors respectively. This is precisely one less than optimal. 
Astley's algorithm consistently gave errors that were smaller than those 
given by Abrahamson's algorithm. The results for the second problem are 
shown in figures 10 and 11. These curves show that Astley's algorithm 
still gave at least the same rates as for the first problem while Abraham- 
son's algorithm converged slower (in the pressure) for this problem, 
demonstrating that Abrahamson's code has convergence rates that are 
problem dependent. 
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2. Non-Uniform Duct Results 

The problems with convergence reported by Abrahamson [7] were 
encountered with the non uniform mean flows that are usually encountered 
with variable geometries. In order to compare the performance of the 
codes in such a situation a quartic half duct was chosen as a test 
problem. The quartic duct under consideration was symmetric about the 
throat and had terminations that were uniform ducts. The throat to 
exit area ratio was chosen as 0.5 and the duct length and termination 
height were both set to 1.0. Three mean flow cases were considered: 

Case (a) The simplest mean flow was constructed by solving the one 

dimensional nozzle equations to yield an x component of mean 
flow velocity, density and speed of sound. The y component 
of velocity at the upper wall was chosen so that the flow 
would be tangential, elsewhere it varied linearly from zero 
on the center line to its maximum value on the upper wall. 

The inlet/outlet Mach, number was chosen to be 0.3 yielding 
a Mach, number of 0.86 on the duct centerline at the throat. 
Case (b) The next mean flow was constructed by modifying the flow 
described in part (a) above to give a boundary layer at 
the throat. To create the boundary layer a weight (w) was 
calculated according to 

w(x,y) - sin [-j(I - (i- x))] 

& 

where 

Y(x) is the duct width at x. 
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The modified velocity field was then the product of the 
velocity field given in part (a) and the weight w. At 
the throat this weight gives a sinusoidal boundary layer. 

The density and speed of sound were not adjusted. 

Case (c) The third mean flow was constructed like the one in (b) above 
except that a much steeper boundary layer was chosen. A 
one-seventh boundary layer approximation, similar to the 
empirical relation observed in pipes, was used. The 
weight (w) chosen to yield such a layer is given by. 


w(x,y) 


[1 


Y(x) ^2 


1/7 

(A- x) ] 


The mean flow, for all three cases, was evaluated at the nodal 
values of the finite elements, and the element interpolation functions 
were used to determine the flow within an element. Eight noded isoparametric 
elements were used for all cases. Case (c) has very steep gradients so the 
element interpolation functions would not approximate these gradients very 
well on a coarse mesh, however, since this is typical of what happens in 
practice we did not make any attempt to correct this problem. The mean 
flow for case (c) has a singularity in its derivatives at the throat so 
the integrals defined for the Astley and Abrahamson weak problems may not 
be continuous. 

Since the mean flow varies rapidly near the throat and at the upper 
wall, meshes were chosen that had a higher density of elements in these 
regions. The meshes were constructed using the following transformation 
for the nodal spacing in the x direction. 
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x(g) - j [i + - g - ( [ + x } " ] > - 1 !?! 1 x - 1.0 

Equal increments of £ were used to yield nodal x values. The following 
trans format ion was used for the nodal spacing across the duct, 

y(n) - y(x) . n [ ^ x + ^ --y — ^3 » o <n£i x = 1.5 

i 

I 

| where 

i Y(x) is the duct width at x. 

Again equal increments of n were used to yield the nodal y values 
at each axial location. Along straight sides the midside nodes were 
always placed in the middle. A typical mesh with 5 elements across the 
duct and 10 elements along the duct (a 5 x 10 grid) is shown in figure 12. 
Three meshes, 5 x 10, 10 x 20 and 20 x 40 elements, were utilized for 
each of the three cases (a-c) described above. The frequency (co) was set 
to 2.0 and the problem of a plane wave incident at x~0 (a^ * 1.0) was 
selected. For this problem the meshes described above would be considered 
very fine since they have many elements per wave length. With realistic 
geometries and frequencies meshes with a similar number of elements per 
wave length would probably exceed the capabilities of most computers. 

As no exact solution exists for the non uniform duct problem a 
comparison of plots of pressure amplitude contours is made. It is expected 
that if the contours are smooth where the mean flow is smooth, and if the 
contour pattern only changes slightly as the mesh is refined , the 
solution is converging . When appropriate, plots of the pressure 
are presented for comparison of the schemes. For the finest grid, 
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calculations were made of the L difference between the pressure given 
by the potential flow code and the pressure given by the other two 
algorithms. 

Figures (13a-c) show the pressure contours for case Ca) given by each 
of the three codes with the 20x40 mesh. Using the criteria discussed above, 
it appears as if the potential flow code and Astley's algorithm are 
converging, however the plot given by Abrahamson's algorithm is not 
particularly smooth and it is difficult to conclude anything. One 
conspicuous difference between the potential flow code and the other 
two algorithms is the apparent rate of convergence (cf. the uniform duct 
solutions). Figures (14a-c) show the pressure contour plots for each 
of the three algorithms given with the coarsest (5 x 10) mesh. 

Figure (14a) shows that the potential flow code contours are qualita- 
tively the same as those for the finer mesh, figure (13a) , however, the 
contours for the other two algorithms show the presence of large wiggles 
which almost obscure any pattern present. 

Figures (15a-c) show the pressure contours given by the three 
algorithms for case (b) on the finest mesh. Case (b) exhibits the 
trends discussed in case (a) above and, if anything, to a greater 
degree. Comparing the contour patterns given by the potential code 
and Astley's algorithm shows a difference in the pressure near the 
upper wall. Figures (16a-c) show the wall pressure plots corresponding 
to the three contour plots given in figures (15a-c) . It is apparent 
that the rotational codes develop a different pressure profile at the 
upper wall near the throat. Another difference between the rotational 
and irrotational codes is the level of the solution. This is indicated 
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qualitatively by the- peak pressures in the throat given in figures 15, 

2 

and quantitatively by the L norms of the solutions given in figure 17. 

Figure 17 contains the norms of the solutions for each of the three 
2 

cases and the L difference between the potential flow pressures and 
the pressures given by the rotational codes for the finest grid. 

Figures (18a-c) show the pressure contour plots given by each of 
the three codes for case (c) . Inspection of these plots shows that 
Abrahamson ' s algorithm almost certainly is not converging and that the 
potential flow solution almost certainly is. The results given by 
Astley's algorithm are inconclusive and are very similar to those of 
Abrahamson' s algorithm for cases (a) and (b) . Figure (19a-c) show 
plots of the pressure along the lines x“0, x=0.5 and x*»I.O given by 
each of the algorithms on the finest mesh. The pressure given by the 
velocity potential code at the throat contains a discontinuity. The 
pressure at the throat given by Astley's algorithm is similar but it 
is reminiscent of the pressures shown in figures (5) and (6) where 
continuous pressures were used to model a pressure jump that was modeled 
much better by the discontinuous elements. The pressures given by 
Abrahamson 's code at each of the three axial locations contain tremendous 
wiggles. 

2 

Figure 17 tabulates the L pressure norms given by each algorithm 

2 

for cases (a-c) on each of the meshes and the L difference between the 
rotational and irrotational pressures with the finest mesh. It is 
apparent that the pressure is consistently larger for the potential flow 
code than the rotational flow pressure and it appears that if the Abrahamson 
and Astley algorithms converge, the pressure is similar. The L difference 


38 


between the pressures given by the irrotational and rotational algorithms 
was greatest for case (b) where they differed by about 30%. The difference 
for cases (a) and (b) were approximately 10% and 15% respectively. 
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40 


Figure 12. Quartic duct profile used in nonuniform duct convergence 
studies. Shown here is the coarsest mesh, 5 x 10. 











PRESSURE AMPLITUDE CONTOURS 



Pressure contours for nonuniform duct case 
Astley-Eversman formulation. 20 x 40 mesh. 
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ONTOUR INCREMENT = O.JOOO 

Pressure contours for nonuni form duct case (c) 
Abrahamson formulation. 20 x 40 mesh. 







Figure 14{a). Pressure contours for nonuniform duct case (a). 

Velocity potential formulation. 5 x 10 mesh. 
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Figure 14(b). Pressure contours for nonuniform duct case (a). 

Astley-Eversnan formulation. 5 x 10 mesh. 




PRESSURE AMPLITUDE CONTOURS 
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Figure 14(c). Pressure contours for nonuniform duct case (a). 
Abrahamson formulation.* 5 x 10 mesh. 


PRESSURE AMPLITUDE CONTOURS 
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Pressure controus for nonuni form duct case (b). 
Velocity potential formulation. ?0 x 40 mesh. 







PRESSURE AMPLITUDE CONTOURS 



Pressure contours for nonuniform duct case (b) 
Astley-Eversman formulation. 20 x 40 mesh. 








PRESSURE AMPLITUDE CONTOURS 
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Pressure contours for nonuniform duct case (b) 
Abrahamson formulation. 20 x 40 mesh. 
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PRESSURE AMPLITUDE CONTOURS 
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Pressure contours for nonuniform duct case (c) 
Velocity potential formulation. 20 x 40 mesh. 
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Pressure contours for nonuni form duct case (c) 
Abrahamson formulation. 20 x 40 mesh. 







PRESSURE P0F1LES ACROSS THE DUCT 



PRESSURE AT X=l + PRESSURE AT 




PRESSURE P0F1LES ACROSS THE DUCT 



58 


PRESSURE AT X = 0 a PRESSURE AT X = 0.5 + PRESSURE AT 
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PRESSURE AT X-0 a PRESSURE AT X = 0.5 + PRESSURE AT 


CONCLUSIONS 


A combination of analysis and numerical experimentation has been 
carried out with the goal of examining the convergence characteristics 
of three finite element algorithms for the modelling of acoustic propaga- 
tion in nonuniform ducts carrying a mean flow. The first model is based 
on the assumption of an irrotational mean flow and irrotational acoustic 
perturbations. This model is currently in use in the calculation of 
propagation in turbo-fan engine inlets and the subsequent radiation to 
the far field. In the present study the physical constraints are relaxed 
to permit a mean flow which is not irrotational, but which is imposed on 
the irrotational model by entering it as data in the model. As viewed 
at this point, this is a non-physical model but serves to broaden the 
scope of the investigation. In the irrotational model the field variables 
are acoustic velocity potential and acoustic pressure Cor density) . 

The other two models do not impose the irrotational assumption. As 
a consequence they use the primitive variables, pressure and two components 
of velocity. One of these models, due to Astley and Eversman [3], takes 
advantage of natural boundary conditions and introduces the idea of 
specifying acoustic input, reflection, and transmission boundary conditions 
in terms of modal amplitudes. The second rotational model is due to 
Abrahamson [6,7] and uses forced boundary conditions. 

In Appendix I it is shown by analysis that under stated restrictions 
the velocity potential formulation will converge optimally. Numerical 
experimentation on both uniform and nonuniform duct models verifies the 
analysis. The stated restrictions are somewhat different than those used 
in current turbo-fan radiation codes in that (1) velocity potential and 
pressure are solved for simultaneously (as opposed to solving for potential 
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and then obtaining pressure by calculation) and (2) the approximation 
for velocity potential is sought in the subspace of continuous functions 
while the approximation for pressure is sought in the subspace of dis- 
continuous functions (as opposed to the subspace of continuous functions 
for both potential and pressure) . 

No analysis is presently available to isolate the convergence 
characteristics of either of the two models based on the primitive 
variables. Information is derived entirely by numerical experimentation. 

It is found that in the examples considered optimal convergence is not 
achieved and in some instances lack of convergence is observed. In 
general, it appears that the Astley-Eversman algorithm displays incon- 
clusive convergence characteristics, certainly not optimal, but not clearly 
nonconvergent . The Abrahamson algorithm is much more likely to display 
nonconvergence and this behavior is clearly problem dependent. 

It was also noted that the potential formulation, with a rotational 
mean flow forced on it, does not converge to the same result as the fully 
rotational models do (when convergence occurs) , the difference being 
more obvious in the boundary layer cases in which a strong shear layer 
exists. Furthermore, the tendency toward nonconvergence in the rotational 
models is exacerbated by a strongly sheared boundary layer. The poorest 
convergence occurred in the numerical experiments using the 1/7 power 
boundary layer. 

It is concluded that: 

(1) The potential flow model converges with reasonable assurance. 

(2) The rotational flow models are subject to slow convergence or 
nonconvergence. 

(.3) Poor convergence or nonconvergence is most likely to occur 
with strongly sheared flows. 
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ABSTRACT 


We investigate a finite element approximation of a problem in 
duct acoustics where the medium is inviscid, the motion is assumed to 
be irrotational and non-homogeneous boundary conditions are imposed 
which arise from the approximation of an infinite duct with one that 
is bounded. A study is made of restrictions for the boundary conditions 
and their approximation for which a finite element approximation of the 
mixed-hybrid type is well-posed. In particular we prove that under a 
reasonable hypothesis on the data that the method will converge 
optimally if the sequence of boundary coefficients has more than h * 
members and the degree of polynomial (k) for the pressure is greater 
than or equal to the degree, k+1, which is used for the velocity potential. 



1 . Introduction 


The problem of interest is the following: Find (<j>,p, {a n >, 

such that 

-2 -2 

(1.1a) -to) c p + V • (pV<J> + c Up) =0 in ft 

(1.1b) p -1 p + U*V<}> + ^iuxf) =* 0 in ft 

(l.lc) (pV4> + c Up)*N = P r 2 <5 (a - a ) cos(mry/d.) 

- - r 1 n=0 n 

(1. Id) -(pV<J> + c" 2 U p)*N » P r 2 5 n^ b n “ b r^ cos ( nn y/d 2 ) 

00 

(l.le) <M r - 2 (a+ + a“) cos (niry/d 1 ) 

1 n«0 

00 

(1 . If) tip - 2 (b* + b”) cos (mry/d 2 ) 

(1,lg) f£ " 0 on r 3 U r 4 

2 4 
ft C 1R represents a duct with boundary 3ft « (j T. where V. ■ {(0,y): y € ]0,d. [} 

i=l 1 

and r 2 * {(L,y): y 6 ]0,d 2 [}. In addition F^ and F^ are assumed to be such 
that ft satisfies the cone condition. The variables of interest p(x), <p(x) 6 C 

oo — 0 

define the acoustic pressure and velocity potential; U € L (ft) represents 

OO — 

the mean velocity; c, p € L (£2) denote the sound speed and gas density; u) 

denotes the pertubation frequency; and {a + } = {a + }* and {b - } = {b“}°° „ 

n n n=0 n n n—O 

represent incoming waves and are assumed to be known. Also, 

6 n = ^ [k F. " d-M? )(mr/d j .) 2 ] Js ; i - 1,2 

where M = | U | /c and k ■ u)/c represent the mach number and characteristic 
wavelength respectively and M„ , etc. denotes its restriction to T., i=l or 2. 
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Equations (l.la-b) are a model for an inviscid, non-heat conducting, 
ideal fluid which is undergoing an irrotational, "small", periodic-pertubation 
about a subsonic mean flow (M < 1) in a duct. The boundary conditions (l.lc-d) 
define compatibility conditions between incident and reflected waves in a 
uniform duct (cf. Eversman and Astley [3]). 

The compressible flow equations have received considerable, attention 
in the literature (e.g. Bristeau et al. [1]) and equations (l.la-b) represent 
a linearization of these. Also Fix and Nicolaides [4] have studied a mixed 
model with Dirichlet boundary conditions for the case when the flow is 
possibly rotational but U » 0 (i.e. a mixed model for Helmholtz's equation). 

In particular they have investigated conditions under which an increased 
convergence rate holds for the pressure (see also Fix et al. [5] and Fix 
and Gunzburger [6 ] ) . 

We will study a finite element approximation to (1.1) which is 
mixed in the sense that (d>»p) are approximated independently and hybrid 
in the sense that (<J>,{a n ) ,{b*}) are approximated independently. In section 2 
we define a weak problem for the boundary-value problem and specify a finite 
element approximation. Then in section 3 we prove that the weak problem is 
well-posed and in section 4 we prove that the finite element approximation 
converges under suitable hypotheses. 
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2. The Finite Element Approximation 

We first define a weak problem which is associated with (1.1) 
and then consider a finite element approximation of this problem. 

First we fix some notation. 

Let H S (fi) , °° > s 0 , denote the usual Sobolev space over 

oo 2 °° 2 2 

the complex field and define S * { {8 n ^ n= Q 6 ^ CMltO}) : Z “ISjJ + |SqI < 00 ) 

n=l 

9 

where -c QNU(0}) denotes the space of square- summable • sequences of complex 
numbers. In addition || • || g ^ and || • ||^ will denote the usual norms for 


H' 


'(ft) and 5 respectively. Also set C with meas > 0, let 


H * {g 6 (Q) : g| r * 0} and define, 

5 


n(n) 


1 n = 0 

n n e IN 


B ^(g n> ^) » ^ | <|>(y) cos (n-ny/d^dy V {g^} 6 S, ip 6 H ; i * 1,2 


n 


Weak Problem : Find (<J>,p»{a m ),{b+}) 6 H x L^(ft) x S x S such that, 

(2.1a) 


-2 - -2 - 

[ - iix. pip + (pV<p + c Up)* Vtf;] 


S2 


+ P r Z <5* (a~i|0 + p r Z 6* B^b+W 

ln*fl "2 n=0 


w uu 

p r Z B*(a^ l>) + p r Z fiV(b“ 4>) 

r i n=0 “ n n r 2 n=0 n n n 


(2.1b) 


[p *p + U.V(j) + ^cio4> ] q * 0 


n 
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(2.1c) 


£ n (n) (a~ + a*)^ - (2/dj) £ i* B*(S n ,?) + dj l Bj(5 Q ,?) 

n=0 n n»l 

oo 00 o . i 2 - 

(2. Id) £ n(n)(b + + bj B„ - (2/d ) £ n B (0 ,$) + d~ B (0 ,*) 

^ n n n l , u u ^ w u 

n=0 n=l 

V 0/;,q,{a n },{3 n }) e H x L 2 (ft) x S x S 

Remark 2.1 . The requirement that <f> * 0 on is a technical detail 

which is dependent on the specific choice of data. For convenience we 
will always include this condition. 

Remark 2.2 . A weak problem with fewer unknown variables may be constructed 
by eliminating the pressure, {a^} and {b*}. However these variables are 
required by Acousticians for comparisons with experiments. 

In constructing the finite element approximation we shall need 
the following standard notions ( cf. Ciarlet [ 2]). Let denote a regular 
family of triangulations of ft where h denotes the maximum diameter of the 
element partition. Define 

ff^(T) * {polynomials of degree < k on T}, T 6 

a /v 

^(T h )» (g 6 C(ft): g|p * 0, g| T » f + i.t where f, f e I» k (T) 

V T 6 T. ), k > 1 

n — 

]M k 1(T h ) = { S 6 l2( ^ ): sl T = f + ^ Where f ’ f 6 ff k (T) 

V T € T t }, k > 0 

n — 

N 

\ - { (g n } N C C: £ n| g n | 2 + |g 0 | 2 < »}, N > 0 

n=0 n=L 
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Here C(ft) denotes the complex-valued functions which are continuous. 
Approximation : Find (<J> h »P h > 6 ^(T^) xM m^V x x 

such that , 


(2.2a) 


r 

[ 

4 


n 


- ‘Lose 2 p, K + (pv*u + C ~ 2 U pj • 


h h 

N 


'h' ; r h j 
N 


+ "r, E „ s 2 •i's'V + "r, J. 4 4 ( KW 

1 n*0 2 n=0 


^11+ N 2 2- 

p_ E St B^Ca; ,*.) + P r 2 r B^(b ,*.) 

r i n -0 ° n n h r 2 n-0 n n n h 


(2.2b) 


(2.2c) 


[p _1 p h + U. V«p + 'U)<J> h ] * 0 


Q 


E S(n)(«: H + a* )S* - (2/d t ) E ■ 

n*0 n-1 


(2. 2d) 


. I i(nHb* H + t; ) 1“ - <j/d 2 ) E n 

n*0 n=l 


V (*i 1 ’V ta n , ’ {a n} > X S N X S B 


- 1 , 


where t > 1, m > 0 and N >_ 0 


2. 1 Preliminaries 


sat || • || 0 = II * |l o n and II • || t - II • || 1>n . we will list some 
lemmas and hypotheses that will be needed in the sequel. First the 
standard Poincare inequality. 

T^mma 2.1. Suppose ft satisfies the cone condition and meas >0. 

Then there exists a X > 0 such that 

Iloilo 1 x II* Ho V* e H 

Lemma 2.2. Let {a } e S, ip 6 H and suppose ft satisfies the cone condition. 
n 

Then for i=l,2 there exists a positive constant yCd^ such that 

Z n(n)|B^(a n ,ijj)| < yd*? max (1 , (d^Tr)* 5 ) || {a n > || s || ip || L 
n=0 

Proof. First note that by the trace theorem. Lions and Magenes [7 ], there 
exists a y > 0 such that V ip 6 H, \#i | € H^r.),. and || |L £ i Y || ip Ik for 

1 i 1 ' s * i i 

i - 1,2 since I\ C 3ft. Then because of the definition of B n (*,0 and 

the Cauchy- Schwarz inequality, B*(a n ,i|j) exists and 

Z n(n)|B^(a n ,<|;)| < dl || (a n }|| s || {<j£}|| s 
n=0 

where ^ = lire*" with {e'*'} being defined by 
n J n n 

r. 



Now by using classical interpolation theory for Sobolev spaces ( cf • Lions 

It 2 Ir 

and Magenes [ 7 ]) * {g: g 6 Dom([I-D ] )} may be constructed as an 

2 1 

intermediate space between L (1*^) and {g 6 H (I\) : Dg(0) * Dg(d^) = 0) • Hence 

|UI|? r * 2 U+(nir/d ± ) 2 ] h |i|£| 2 > min(l,jr/d ) || {*£} |£ 
x n=0 x n x n i 

which leads to the desired estimate after we combine inequalities. 

Next we cite a standard result of approximation theory ( cf. 

Ciarlet [ 2]). Let I • I n , for m 6 IN, denote the usual seminorm on H 01 (ft) . 

Lemma 2.3 . Suppose ip 6 (ft) , k ^ 1 , and q € H m+ * (ft) , m _> 0 . Then if 
is a regular family of triangulations of ft, there exists a Cp Cj > 0 
(independent of h, t|> and q) such that 


inf II * - v h Hi - c i hk H’lk+i.n 


v h e W 


i“f ll <1 - v. ||. < C,h m+1 |ll 






We now list some reasonable hypotheses that the known variables 
should satisfy. 

Hypothesis 2.1 . {a*}, {b^} 6 S ; c, p € L (ft) ; w 6 {g 6 IR: g 0}; 

and U € L°° (ft) 2 with 

ess inf c(x) = c Q >0 
x 6 S 
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II C Hl”«!) ’ C 1 


ess iii£p(x) * Pq > 0 

X 6 3 


II p llwa) " P 1 


II U II _ 2 = u 

- L“(fl) 


N 

Here L “(ft) , for N ■ 1,2 is defined on the real scalars. 

We close with a standard approximation result for sequences. 
First we define for A > 0 


S* * {{g n ) 6 £^QN U {0}): Z n(n) 

n=0 


2a. + 1 1 1 2 . i 

|gj < -} 


Lemma 2.4. Let {a } 6 A > 0. Then for N > 0 there exists a positive 
n — — 

scalar C such that 


inf || {a n - g n }|| s < C N -/l |{a n }|| 

{g } 6 S 

n N 
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3. Existence and Uniqueness of Solutions 


2 2 

We note that the form on (H x L (fl) x S x S) , which represents the 

left hand side of (2.1a-b) and the right hand side of (2.1c-d), is sesquilinear 

2 

and the form on H x L (ft) x S x S, which represents the right hand side of 
(2.1a-b) and the left hand side of (2.1c-d), is antilinear. Thus to prove 
that the weak problem has a unique solution it suffices to show that the 
conditions of the Lax-Milgram theorem are satisfied. 

First observe that from the definition of 6^, along with Hypothesis 

tl 

2.1, that 

(3.1) |5*| < fe ± n(n) for i - 1,2 

where k i * [k 2 + | 1 - M 2 | (ir/d^) 2 1" 5 ; let u » (<j>,p,{a n ), {b*}) » 

i i 

v - (<|>,q,{a n M8 n }), U - H x L 2 (ft) x S x S, || u ||Jj * || <(> || 2 + || p || 2 


lK a n }|ls + il^n^^S Vu 6 U; and define 


+■. ii 2 


F(v) - p r Z Sj- B* (a+ip) + p f Z 6 2 B 2 (b“((i) 
r i n«0 n n n r 2 n-0 n n n 


- Z n(n) (a* a + b~ £ ) 
n=0 n n n n 


A(u,v) 


-2 - -2 _ - I - 

[ - -Uii c pip + (pV<p + c U p) • Vip + p pq 


ft 


1 T, 


+ (U*V(j))q + -tuxjiq - (2/d.) Z n B (a. ,<j>) + £ n(n)a a 

1 n-1 n n n=0 n n 


- d“ l Bj(oL,5) - (2/d 2 ) Z n B* (B ,5) + Z n(n)b* £ n 

n=l n=0 


- d 2 lB o ( ' 6 o-* ) + B tPr«i + °r { n B »h» 

n=0 1 2 
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Lemma 3.1 . Suppose Hypothesis 2.1 holds and SI satisfies the cone 


condition. Then F 6 II' and A is a continuous sesquilinear form on U x U. 

Proof . From the definition of F, use of the Cauchy-Schwarz inequality. 
Hypothesis 2.1, 3.1 , and Lemma 2.2 we have that 

|F(v)| £ P r ^ fejYd^ max (1 , (dj/ir)' 5 ) || {a+} |j s || iff || L 

+ p r 2 fe 2 yd 2 max || {b”) ||^ Ikllj 

+ IK<> Hs H { «n } Hs + ll {6 n } Hs 

thus 

II F ll u - i max/l, max[p r d^fe i max (1, (dj/iO^yA ( || (a*} || 5 + ||{b~>|| 5 ) 

' i” 1)2 ' 

hence F 6 U ~ . Now if we perform similar operations on A(*,*) 

/ —2 -2 —1 

|A(u,v)| £ maxi me” , Pj, UjCq , p” , Up <u, l, 

max [ (Pj. k.jd*? + Zd^* 5 ) max(l , (d^/ir)" 5 ) ]Y^|I u|l g|l v II (j 
i-1,2 1 / 

which completes the proof. 

Define ^ » y max [ (2 + ^d^d”* 5 max (1 , (d^ir)" 1 ) ] and + w//A) (l+c Q 2 ) • 

i-1,2 i 

Theorem 3.1 . Suppose Hypothesis 2.1 holds and SI satisfies the cone condition. 
Then if there exists a positive number e such that 

2 < 8 e <C4P 0 A/(1+A) - 4 P 1 )/5 1 

there exists a unique u e U satisfying (2.1) with || u ||y £ C (||{a^}||^ 

+ ||{b~) ||^) 'where C > 0. 
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Proof . Lemma 3.1 holds therefore we need only prove that A is coercive 


on U x U. Using Lemmas 2.1 and 2.2, Hypothesis 2.1 and the Cauchy-Schwarz 
inequality. 

Re A(u,u) >_ o 0 X/(l + A) || * ||j + p, 1 1| p ||jj + ||fc n I||| + ll<>>*>ll| 

- C 2 II P Bo II * 111 -C l ll*lli(K , lls + ll {b a ) lls) 

and by applying the arithmetic-geometric mean inequality to the last two 
terms we have with * PqA/(1 + A) - - 25^, C 2 = 1/Pj - 

and * 1 - 5^/4e: 

Re A(u,u) > C L || 4> || 2 + C 2 || p ||q + C 3 (||{a”>||| + ||{b+}|| 2 ) 

Now e > 0 may be chosen so that 

C 1 "" C 2 = P 0 X/(1 + X) ' [C 4 + (C 4 + ^2 )1Sl/2 _ 25 l e 

where = pgA/(l + A) - 1/p^ - 2£^e. Thus by hypothesis (after some 
manipulation) C^, > 0. Hence A is coercive on U x U and the Lax-Milgram 

theorem holds therefore completing the proof. 

Theorem 3.1 may also be extended to apply to the approximation. 

First note that W^CT^) with || • ||j and (T^) with || • ||q are closed sub- 
spaces of H and L 2 (fl). Furthermore if we extend {g|J} 6 to {g^} 6 S by 
defining = 0 V n > N then with || • ||^ is a closed subspace of S. 

Thus by repeating the proof for the finite dimensional case we have: 

Theorem 3.2 . Suppose the hypothesis of Theorem 3.1 holds. Then there exist 

a unique u^ = ^h’ P h* * a n N * ’ ^ b n N ^ 6 U h = x]M ra 1 ^ T h^ X 5 N X S N satis “ 

fying (2.2) with || u h || u <C (||{a+ }|| s + ||{b~ }|| s ). 
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. Convergence of the Approximation 


We now prove convergence of the approximation (2.2) to the solution 
of (2.1) under an hypothesis that the solution and the data is "regular". 
Set v h - (i|; h ,q h »{a^},{3^}) 6 and define 


4 „1, + 


2 „ 2 , - 


’ (v ) - p Z 6 B (a , ^».) + P r Z 6 B (b ,ip.) 
N h F 1 n-0 n n “ h 2 n«0 n n n h 


4- —NT - -M 

- Z n(n) (a a + b B ) 
n-0 n n n n 


A N (u h ,v h ) » | [ - -auc” 2 p h t(l h + (pV<|> h + c 2 U p). • V\p + p L p h q t 


1 /-N T 


-N -N 


+ (U*V<j> h )q h + iw<|) h q h ] - (2/dp Z n B n ( 0 n *+ h ) + s u ( n ) a n a n 
~ ~ n“l n=0 


- d i‘ B i <5 0-*h ) - (2/d 2> ", " B a + '‘< n > b n N ' 8 " 

n*l n*0 

,-1^2 t \ ~ P f l _1, - N . N ~2 _2 + N . 

- d 2 WV + £, tP r , B n (a n 'V + B r. (b n ■ V h )1 

n-0 1 2 

Theorem 4.1 . Suppose the hypothesis of Theorem 3.1 holds, is a regular 
family of triangulations, £ > 1, n > 0, A > 0, u and u^ are solutions to (2.1) 
and (2.2) with u 6 (Hfl H^ +1 (fl)) x H° rt ' 1 (fl) x S* x S'* and { a *Mb n > 6 ^ * 


Then there exists a positive constant C such that. 


“ - «h fc + ^‘iPi^X.a * H '' l( H U n 1 IU 


+ « { s>ii^ + K>i^ + iit b ;>i^) 
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Proof . Let e * u - u h * Then A(e,e) - A(e,u-v h > + ACe.v^u^) \J v h 6 
since A is sesquilinear. But A is continuous and coercive on U x U 
(and x (i^) ; u and u^ are solutions to (2.1) and (2.2); and A 1J (u^,v h ) * 
A(u^,v^). Thus 

ii » - “ h iiu i c ni - - \ iiu + i F <vv - F u (u h ■ vi u ■ u h iy 


Now 


F( V " W 


Z [p_ «i“Bj(a^ f *. ) + P r B^(b" * )] 
n-IH-1 r i nnnh r 2 nnnh 


Therefore by using (3.1), noting that V {g Q } e S 
and i » 1 , 2 


Z n(n)|B*(g )| <.N 

n=»N+l n n h 


K Z n(n) |B^(q_(n)' rL g n ,i^ h ) 
n-0 


and applying Lemma 2.2 we get that 



We obtain the desired result after choosing v 'so that Lemmas 2.3 and 2.4 

ft 

hold. 

As might be expected we may not approximate (2.1) with arbitrary 
values for t, m and N without possibly sacrificing accuracy. A judicious 
choice may be realized by observation of the estimate || u - ||y. 

Corollary 4.1 . Suppose that m >_ t~l and N ^ h . Then 

II “ - “h llu * °«>*> 
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Remark 4.1. Corollary 4.1 allows one to determine how many boundary 


terms must be included in the approximation depending on the regularity 
of the boundary data. In many acoustics problems we find that N = 1 
implying that N > (1/h). 


-A- 14- 


4. 1 Numerical Example 

We consider a problem which is not explicitly covered by our theory 
but nevertheless can be analyzed in a similar manner. Let ft = {(x,y): 
0<x<l,0<y< 0.5}, c(x,y) = 1. \J (x,y) 6 ft and define for 0 <_ y <_ 0.5 


P(x,y) 


1 - 


1. 0 <_ x < 0.5 

0.5 0.5 < x < 1. 


Also let e define a unit vector in the x-direction, set 
~x 


U 




-0.25 e 0<x<0.5 
-x — 


-0.5 e 0.5 < x < 1 . 
~x — — 


co * 1 . , a+ * 1. , a + ■ 0. V n > 1 , b” ■ 0 V n > 0 and append to (1.1) 
On — n — 


the following jump condition for x ■ 0.5 ,0<y<0.5. 


( 


lim ( [p(x-e,y)V<Kx-e,y) + c“ 2 p(x-e,y)U(x-e,y) ] 


■u) 


e -*• 0 


( 


-2 


* lim ( [p(jH-e,y)V<Kx+e,y) + c p(xfe,y)U(x+e,y) ] 
£ -*■ 0 


•-) 


We then have a classical problem which has as a solution for <p, {a^} and {b^} 


<j>(x,y) 


exp (- 4/3 -6c) + a n exp (4/5 6x), 0 < x < 0.5 


< 


exp (- 2 6[x-l]), 0.5 < x < 1. 


a n = 0. 161 - 0.219 -t ; a =0 V n 6 IN 
0 n 


bj = -0.1276 - 1.3272 -t ; b* = 0 V n 6 K 


The pressure may be calculated using (l.lb) and the jump condition. 
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Remark 4.3. This solution does not satisfy the hypothesis of Theorem 4.1. 


However the same convergence rates will hold for our approximation as long 

as the graph {(0.5,y): 0 y 0.5} does not intersect the interior of any 

element in T, . Otherwise, as is well known, we would have to use discontinuous 

h 

interpolating functions for the approximation. 

Three n ume rical experiments were performed using uniform rectangular 
elements. We denote a uni form grid with I elements in the x direction and 
J elements in the y direction as a "J x I grid". 

The three cases we consider are: 

(1) (\,Ph) 6 x with 5 x 10, 10 x 20 and 20 x 40 grids. 

(2) (4>k»Ph) S ^ x IM^with 5 x 11, 10 x 21 and 20 x 41 grids. 

(3) ^h’^h^ 6 x 5 x 10, 10 x 20 and 20 x 40 grids. 

The error in pressure is tabulated in Table 1 for each of the 
experiments. Experiment 1 yielded the desired convergence rate (rate « 2) 
while the other two experiments did not (rate ~ 1) as expected. 


-A-16- 


References 


[1] M. 0. Bristeau, R. Glowinski, J. Periaux, P. Perrier, 0. Pironneau 

and G. Poirier, Transonic flow simulations by finite elements and 
least square methods. Proceedings of the Third International 
Conference on Finite Elements in Flow Problems, Banff, Alberta, 
Canada, 1980, pp. 11-29. 

[2] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, 

North-Holland, Amsterdam, 1978. 

[3] W. Eversman and R. J. Astley, Acoustic transmission in non-uniform 

ducts with mean flow; Part 1: The method of weighted residuals, 

J. Sound Vib., 74 (1981), pp. 89-101. 

[4] G. J. Fix and R. A. Nicolaides, An analysis of mixed finite element 

approximations for periodic acoustic wave propagation, SIAM J. 
Numer. Anal., 17 (1980), pp. 779-786. 

[5] G. J. Fix, M. D. Gunzburger and R. A. Nicolaides, On finite element 

methods of the least squares type. Comp. & Maths, with Appls., 

5 (1979), pp. 87-98. 

[6] G. J. Fix and M. D. Gunzburger, On numerical methods for acoustic 

problems. Comp. & Maths, with Appls., 6 (1980), pp. 265-278. 

[7] J. L. Lions and E. Magenes, Non-Homo geneous Boundary-Value Problems 

and Applications, Springer, Berlin, 1972. 



























